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SUMMARY 


A ccxnputer program for the calculation of the power-density spectrum (PDS) 
fran a data base generated by an Advanced Continuous Simulation Language (ACSL) 
simulation program is presented. The program uses an algorithm that employs the 
fast Fourier transform (FFT) to calculate the PDS of the variable. This is done 
by first estimating the autooovariance function of the variable and then taking 
the FFT of the smoothed autocovariance function to obtain the PDS. The PDS pro- 
gram is written so that it is transparent to the ACSL run-time executive (that 
is, the time history of the recorded variable is replaced with its PDS) and 
so that the ACSL user can perform a frequency analysis from the time-domain ACSL 
simulation model. An example of the use of the program is presented that deter- 
mines the frequency content of the output waveform for a Van der Pol oscillator. 


INTRODUCTION 

In engineering science, mathematical models are often developed to aid in 
the understanding of the dynamic behavior of the physical system. These models 
also are of benefit in examining the results of proposed modifications in the 
physical system. The engineering models are often developed with the aid of 
digital computer simulation programs. One of these programs is the Advanced 
Continuous Simulation Language (ACSL) described in references 1 and 2. The ACSL 
system is designed for modelling the behavior of continuous systems described 
by time- dependent, nonlinear differential equations and transfer functions. The 
inputs into ACSL are in two distinct groups: the first group contains those 

inputs concerned with defining the model or structure of the system being simu- 
lated; the second group contains the sequence of commands that execute this 
model - i. e. , change parameters, start runs, plot data. One of the outputs of 
the ACSL program can be the time history of specified dynamic variables, sampled 
at regular intervals and stored on a mass storage device. This record may be 
compared with the behavior of the physical system to determine the fidelity of 
the mathematical engineering model. 

The accuracy of the engineering model may also be verified by use of power- 
density spectrum (PDS) analysis. The PDS of the model can be calculated and 
then compared with an experimentally measured PDS of the physical system. One 
method used in the past to do this analysis utilizes a numerical quadrature pro- 
cedure, to calculate the complex Fourier integral, at a number of selected 
frequencies (see example program number 9 in ref. 1) . If a complex system were 
being analyzed or a large number of frequencies were required to define the PDS, 
the numerical quadrature method could use a substantial amount of computer 
resources. This paper describes a computer program which utilizes a fast Fourier 
transform (FFT) technique to estimate the PDS of the model. Because of the 
efficiency of the FFT, this program should require less computer resources than 
other methods. The program is designed to be an extension of the ACSL system so 
that the PDS analysis can be performed as part of the ACSL sequence of commands 
that exercise the model. In the remainder of this paper, the PDS algorithm used 



in this computer program is described, the computer program's subroutines are 
discussed, the use of the program is described, and an example of the program 
execution is presented. 


PDS ALGORITHM 

In this section, an overview of the algorithm used to calculate the PDS is 
presented. This section is not intended to be a tutorial on the subject of PDS 
analysis and computation. For further information on this subject the reader 
should consult reference 3. 

The PDS S(0]!) of the continuous dynamic variable x(t) is defined as the 
Fourier transform of its autocorrelation function R(t), or 


S(co) 



R(T) dx 


where the autocorrelation is defined as 


R(T) = E[x(t)x(t+T) ] 

and E is the statistical expectation operator. Note that if x(t) is a zero- 
mean function, i.e., E[x(t) ] = 0, then R(x) is also the autooo variance 
function. If a sequence of numbers is derived from the continuous variable by 
periodic sampling, the Fourier transform of the continuous variable can be 
approximated by 
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where S(fi) is the discrete Fourier transform (DFT) of the sampled-data autocor- 
relation function and T is the sampling period. The FFT algorithm is a fast 
method of computing the DFT at a discrete number of frequencies. The FFT is used 
by the program discussed in this paper not only to calculate the PDS from the 
autooovariance sequence, but also to calculate the autooovar iance sequence frcm 
the input data sequence. 

The PDS algorithm used by the computer program described in this paper is 
suggested by Oppenheim and Schafer in reference 3. An estimate of the PDS is 
found taking the FFT of the smoothed truncated autocovariance function of the 
time history. The truncated autooovariance function is calculated by taking the 
time history, composed of a data sequence of N data samples, and dividing it 
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into K data sequences, each of length L/2. The FFT of each of these smaller 
data sequences, augmented with L/2 zero samples (forming a data sequence of 
length L) , are then used in a formula that uses the property of circular con- 
volution to force the values of the autocovariance function to be correct for 
the first L/2 data elenents. In doing so, the truncated autocovariance func- 
tion is calculated by using K FPT's of length L and one inverse FFT for a 
real symmetric data sequence. This autocovariance function is then smoothed 
with a data window and the PDS is calculated by taking the FTT of the real sym- 
metric data sequence. By using this method, N can be selected as large as 
necessary to obtain a good estimate of the autocovariance function, while the 
values of L and the sampling period are chosen to determine the frequency 
range and resolution of the PDS. 

The subroutines that are used to calculate the FFT were obtained from refer- 
ence 4, sections 1.2 and 1.3. In order to gain maximum efficiency, special pur- 
pose FFT algorithms were used for real input sequences and for real, symmetric 
input sequences. The basic FFT routine, which is common to all the special pur- 
pose routines, is a fast radix 8-4-2 algorithm which performs as many base-8 
iterations as possible and then performs one base-4 or base-2 iteration, if nec- 
essary. Additional information on the FFT subroutines used can be found in 
reference 4. 


PROGRAM DESCRIPTIOl 

ACSL is an application oriented system for investigating the dynamic behav- 
ior of jiiysical systems described by sets of differential equations. The inter- 
face of the PDS program to ACSL is illustrated in figure 1. Initially the model 
description is read by the ACSL translator. The output of the translator 
is passed to the run-time executive, which also reads the run-time commands from 
the user. With the proper commands, the ACSL executive generates a data base, 
the PREPARE file, consisting of the time histories of the specified dynamic 
variables. The ACSL executive uses this PREPARE file to generate listings and 
plots of the simulation. The PDS program of this paper, which is called by the 
executive, reads the PREPARE file, calculates the PDS's of the variables, 
and writes them on the PREPARE file. After the program has returned to the ACSL 
executive, the executive can be directed to list or plot the PDS's. 

The structure of the PDS program is given in figure 2. There are four sub- 
routines that were written for this report (PDSMAIN, PDS, WINDOW, INOUT) . A 
listing of each of the routines described in this section can be found in appen- 
dix A. These subroutines call the FFT subroutines FFA, FFTSYM, and IFTSXM, 
which are described in reference 4. The subroutines ZZLINE and ZZPAGE are ACSL 
subroutines for output formating. The subroutines described in the subsequent 
sections make use of the imaginary part of a declared complex number being 
stored sequentially in memory after the real part of the complex number. For 
this reason, care should be taken with program modifications not to disturb the 
storage allocations in the calling statanent to subroutine PDS. In the sections 
which follow, each of the programs developed in this report are discussed. 
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Subroutine PDSMAIN 


Usage ; 

CALL PDSMAIN(A,MP,DELTT,L02,IWIND,IF02,AMAX,IDIV,L0G) 


Argument list ; 

A Work array with dimensions of AMAX or greater. 


MP Equal to the number of dynamic variables written on the PREPARE file 

by the run-time executive. 

DELTT The sampling interval at which the data were recorded. 

L02 The number of data points in the calculated estimate of the autooo- 

variance sequence, L/2. Must be a power of 2. 


IWIND Selects the type of data window to be applied to the autooovariance 
function; 


IF02 

AMAX 


IWIND 

Window type 

1 

Rectangular, or no window 

2 

Bartlett 

3 

Hanning 

4 

Hamming 

5 

Blackman 


The number of data points in the calculated PDS. 

The length of the work array A. Must be greater than; 


512 + 2(MP-1) + MAXIMUM(A,B) 


where 

A = (L02+2) (MP-1) 

B = (IF02+1) (MP-l)+IF02+2 

IDIV If N is not equal to zero, the Nth dynamic variable is divided 

into the remaining variables processed by the program. If N is 
equal to zero, nothing is done. 
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LOG 


A logical variable, that if set true; the logarithm of the frequency 
is placed on the PREPARE file. If LOG is true and IDIV is not equal 
to zero, then the PDS of each variable, except for the Nth vari- 
able, is converted to decibels. 

This routine is called the ACSL simulation program. The purpose of this 
routine is to organize the data storage required by the routine PDS. If the 
amount of storage required exceeds the value of the variable AMAX, subroutine 
PDS will not be called. Instead, the storage requirements will be printed and 
the program will be abnormally terminated. If there is enough storage in 
array A, subroutine PDS will be called. On return fron the routine PDS, the 
amount of storage required is printed along with the mean and variance for each 
of the variables processed by subroutine PDS. PDSMAIN will then return normally 
to the calling program. 


Subroutine PDS 

Usage ; 

CALL PDS(BUFF,MP,KBUFP,A,X1,X1P1 ,MPM1 ,DELTT, XM,X2M,FREK,LP2,IWX,L02P1 , IF02P1 , 
IDIV,RA,RX1,RX1P1 ,FREKX,LOG) 

Argument list ; 

BUFF Data storage buffer used by INOUT with dimensions MP by KBUFF. 

MP Number of variables on the PREPARE file. 

KBUFF Dimension of BUFF. 

A A ccmplex work matrix with dimension MPMl . 

XT A complex work matrix with dimensions L02P1 by MPMl. 

XI Pi A complex work matrix with dimensions L02P1 by MPMl. 

MPMl MP-1. 

DELTT Sampling interval. 

XM A vector of length, MPMl, that on return contains the mean of the 

dynamic variables. 

X2M On return contains the variance of the dynamic variables. 

FREK A work matrix with dimensions IF02P1 by MPMl. 

LP2 L4-2. 

IWX The window-select integer, same as IWIND in subroutine PDSMAIN. 

L02P1 L/2+1 . 
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IP02P1 IF02+1, where IP02 is defined by PDSMAIN. 

IDIV See subroutine PDSM2VIN. 

Kh Work matrix that must be set equal to array A in the calling program. 

On returning to the calling programf the scaled auto covariance 
sequence is returned. 

RX1 A work matrix that must be set equal to the array XI in the calling 

program. 

RX1P1 A work matrix that must be set equal to the array X1P1 in the calling 
program. 

FREKX A work vector with dimension of IF02P1. 

LOG See subroutine PDSMAIN. 


Subroutine PDS calculates the algorithm proposed on page 560 of reference 3. 
The major steps of the routine are as follows: 

1 . The mean and variance of the variables that are recorded on the PREPARE 
file are determined and stored. 

2. The number of data sequences of length L02 that exist on the PREPARE 
file are determined and stored in the variable R. The PREPARE file is then 
rewound. 

3. For each of the K segments, as the data sequence is inputed the mean 
is subtracted and the data sequence is augmented with L02 zero-data samples. 

4. The FFT for each of the segments is calculated by the routine FFA and 
conbined to form the input to the routine IFTSYM, which calculates the estimate 
of the autooovariance function. 

5. The autooovariance function is normalized by the variance and the 
selected data window is applied by the routine WINDOW. 

6. The scaled PDS is calculated, by the routine FFTSYM, by taking the FFT 
of the smoothed autooovariance function. 

7. The PDS is corrected for scaling and outputed to the PREPARE file. 


Subroutine WINDOW 


Usage : 

CALL WINDOW (A,LP2 ,MPM1 , IWX,FREK, IF02 PI ) 
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Argument list ; 


A The input autocovariance sequence that is to be attenuated. 

LP2 The number of data points in the autooovariance sequence. 

MPM1 Number of dynamic variables being transformed. 

IWX The type of vdndow to be used. See subroutine PDSMAIN. 

FBER Output smoothed autooovariance sequence. 

IP02P1 The number of data points in the output sequence. 


The purpose of this subroutine is to attenuate the autooovariance sequence 
with a data window. On returning to the calling program, the array FREK con- 
tains the attenuated sequence of the first L02P1 data elenents with the remaining 
IF02PI - L02P1 data elements equal to zero. 


Subroutine INOOT 


Usage : 

CALL INOUT (BUFF ,N ,M , IN , UN , lERR) 
Argument list: 


BUFF 


N,M 


IN 

UN 

lERR 


Input data buffer. 

Dimensions of BUFF, N*M must be less than 512. 

If equal to 1, the program will read data frcm the input device; if 
not equal to 1, the program will write to the output device. 

Input/Output device number. 

Positive if I/O error occurred; zero if an end-of-file condition 
occurred; and negative if no I/O error occurred. 


Subroutine INOUT is used to transmit data to and frcm the ACSL run-time 
executive. The routine INOUT is used to read and write data to the PREPARE 
file (unit number 8) . The routine will be alxiormally terminated if more than 
512 words are to be transmitted or if an input/output error occurs. If, at 
seme time in the future, the run-time executive's method of writing data on the 
PREPARE file is changed, this subroutine must be modified to conform to the new 
accession method. 
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PROGRAM OPTIONS AND USAGE 


The purpose of this section is to describe several of the options that are 
available and to give some general guidelines on the use of the program. These 
guidelines are as follows: 

1. Since the PDS subroutine replaces the first variable that is stored on 
the PREPARE file with the frequency^ the independent variable of the simulation 
(time/ T) , should be placed as the first variable on the file. 

2. The sampling interval is the communication interval (CINT) in the ACSL 
program. The sampling interval should be small enough so that aliasing does not 
occur. A general rule is to choose a sampling interval at least equal to the 
period of a frequency that is 20 to 40 times the highest frequency of interest. 
This rule is arbitrary because it assumes that the actual PDS is almost zero at 
this frequency. 

3. The length of the data base is chosen so that the estimate of the auto- 
covariance is ergodic. Since this length is dependent on the systen and the 
randan process, it must be chosen by trial and error or experience. 

4. The frequency resolution of the PDS is determined by the sampling inter- 
val (DELTT in the listing) and the number of data points in the autooovariance 
function estimate (L02) and is equal to 1./ (2.*DELTT*L02) . The variable L02 
must be a power of 2 , but since the FPT routine is a radix-8 type, L02 should 
be a power of 4 for the optimum use of the PPT routines. 

5. The PDS may be calculated at closely spaced frequencies, as is conve- 
nient and practical, by using the variable IF02. The frequency increment between 
PDS data points is equal to 1 ./(2.*DELTT*IF02) . A general rule is to make IF02 
at least twice L02. 

6. No matter how large the variable IF02 is, the frequency resolution of 
the PDS is still determined by the length of the autooovariance estimate 
(DELTT*L02) . If in the PDS program, IF02 is smaller than L02, L02 is used 
instead of IF02 for the frequency output interval calculation. 

7. In the development of the program, the only feasible means for input- 
ting the number of variables that are written on the PREPARE file (MP) to the 
program was as a passed variable in the calling sequence. If the user does not 
insure that MP is correct, unpredicted results can occur. 

8. In parameter identification, it is useful to conpare the magnitude of a 
frequency response curve to the PDS of one variable divided by the PDS of another 
variable. This can be done in the PDS program by setting the variable IDIV equal 
to N, where N represents the PDS of the variable that is to be divided into 
the other PDS's. No division takes place if IDIV is zero. 

9. If the logical variable LOG in the calling sequence is set true, the 
logarithm of the frequency is written onto the PREPARE file. Also, if IDIV is 
not equal to zero, the PDS's of variables, except for the PDS of the Nth vari- 
able, are converted into decibels. 
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10. So that the calculated PDS can be easily compared to the magnitude of a 
transfer function, the square root of the PDS is taken before it is written on 
the PREPARE file. 


EXAMPLE OF PROGRAM USE 

An example is discussed in this section which shows the subroutine being used 
to calculate a PDS of a time history and the ACSL run-time executive plotting 
the results. The example presented is an ACSL program simulating a Van der Pol 
oscillator. The equation for the Van der Pol oscillator is similar to that of 
a free vibration of a spring mass system with viscous damping. However, the 
damping term of this equation is nonlinear in that it depends on the amplitude 
of the oscillation. At small amplitude levels, the damping is negative and the 
amplitude of the oscillation grows. At large amplitude values, the damping is 
positive and the amplitude diminishes. The combination of these two effects 
leads to a stable limit cycle. The equation for describing a Van der Pol oscil- 
lator can be written in parametric form as 


X - X (1 - x^)x + X = 0 


where X determines the growth rate of the oscillation and x is the state 
variable that represents the output of the oscillator. A dot over a symbol indi- 
cates the derivative with respect to time. In reference 5, Rogers and Connolly 
give an approximate solution for the fundamental frequency cJq in radians per 
second as 


X2 



The amplitude of the fundamental frequency and the first two harmonics are 
defined as 


Ao = 2 
A-| = 0 

2X 


A listing of the ACSL translator input for this problem is given in appendix B. 
Except for the additional statements needed to specify the calculation of the 
PDS, the simulation program is identical to the example presented in reference 2, 
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The initial cx>nditions are such that with X » 1, the simulation will oscillate 
in a stable limit cycle with a fundamental frequency of 0.14 Hz. The subrou- 
tine is configured to calculate the autocovariance function with data sequences 
of 512 data samples, uses a Blackman data window and calculates a PDS having 
1024 data points. Also, the program expects to find three variables recorded on 
the PREPARE file every 0.063 sec. Since this results in the recording of approx- 
imately 800 data points on the PREPARE file, there will be only one data sequence 
used in the algorithm. For this case, this will be sufficient since the autoco- 
variance is a periodic and deterministic function. With the preceding values, 
the frequency resolution of the PDS will be 0.016 Hz, and the frequency increment 
between the plotted data points of the PDS will be 0.008 Hz. 

A listing of the ACSL run-time executive input can be found in appendix C. 

The simulation is constructed so that the initial START command generates the 
data base to be used by PDSMAIN. This data base is displayed by the next two 
plot commands. Figure 3 is the output generated by the ACSL run-time executive 
for the first plot command listed in appendix C. This instruction produces a 
plot of the variable X, the oscillator output, on the y-axis against the vari- 
able T, time. A phase- plane plot of the system can be generated by plotting 
the derivative of x (XD in the simulation) , against the variable X. Fig- 
ure 4 was generated by the second PLOT command in appendix C. When the second 
START command is issued, the simulation has been reconfigured to call PDSMAIN 
and return to the executive. The last PLOT instruction produces figure 5, 
the desired PDS of the output variable, PDSX plotted against frequency W in 
Hertz. The peak of the first harmonic is not exactly what was predicted in 
reference 4. The applied window reduced the low-frequency peak. However, the 
third harmonic peak is in agreement with reference 5. 


CONCLUDING REMARKS 

A computer program for approximating the PDS from time histories generated and 
stored by an advanced continuous simulation language (ACSL) simulation program 
has been presented. The method uses an algorithm that employs the fast Fourier 
transform (FFT) of the time history to calculate an estimate of the autocovari- 
ance function of the time history. The power-density spectrum (PDS) is then 
calculated by taking the FFT of the windowed autocovariance function. An exam- 
ple was presented that determined the PDS for a Van der Pol oscillator. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
May 15, 1981 
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APPENDIX A 


PROGRAM LISTING 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE PDSMAIN(A,MP ,DELTT,L02 , IWIND , IF02 ,AMAX, IDIV ,LOG) 

PURPOSE TO CALCULATE THE PDS OF VARIABLES 
WRITTEN TO THE PREPARE FILE 


ANALYSIS H.J. DUNN - NASA, L ARC 
WRITTEN MARCH, 1980 

ARGUMENT LIST 
A WORK ARRAY 

MP NUMBER OF VARIABLES WRITTEN ON THE PREPARE FILE 

DELTT SAMPLING INTERVAL 

L02 THE NUMBER OF DATA POINTS OF THE ESTIMATED 

AUTOCOVARIANCE FUNCTION (MUST BE A POWER OF 2) 

IWIND SELECTS THE DATA WINDOW TYPE, 

IWIND WINDOW 

1 RECTANGULAR 

2 BARTLETT 

3 HANNING 

4 HAMMING 

5 BLACKMAN 

IF02 NUMBER OF DATA POINTS OF THE PDS (MUST BE A POWER 
OF 2) 

AMAX LENGTH OF THE WORK ARRAY A 

IDIV IF IDIV=N THEN DIVIDE THE PDS OF THE NTH VARIABLE 

INTO THE PDS OF THE REMAINING VARIABLES ON THE PREPARE 
FILE 

LOG IF LOG IS SET TRUE IN THE CALLING PROGRAM, THE 

LOGARITHM OF THE FREQUENCY IS PLACED ON THE PREPARE 
FILE 

IF IDIV.NE.O AND LOG IS TRUE THE REMAINING VARIABLES ARE CHANGED 
INTO DECIBELS 

DIMENSION A(AMAX) 

INTEGER AMAX 

LOGICAL LOG 

MPM1=MP-1 

L=L02*2 

L02P1=L02+1 

IF02P1=IF02+1 


n 
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IF(L02P1.GT.IF02P1) IF02P1 = L02P1 

LP2=L+2 

KBUFF=512/MP 

Nl=l 

N2=N1+KBUFF*MP 

N3=N2+MPM1 

N4=N3+MPM1 

N5=N4+LP2*MPM1 

N61=N5+LP2*MPM1 

N71=N61+LP2*MPM1 

N62=N5+IF02P 1*MPM1 

N72=N62+IF02P1+1 

N7=AMAX0(N71,N72) 

IF(N7.GT.AMAX) GO TO 100 

CALL PDS(A(N1) ,MP,KBUFF,A(N4) ,A(N5) ,A(N61) ,MPM1,DELTT, 

1 A( N2 ) , A( N3 ) , A( N5 ) , LP2 , IWIND , L02P1 , IF02P 1 , IDIV , 

2 A(N4) ,A(N5) ,A(N61) ,A(N62) ,LOG) 

100 CALL ZZPAGE(1,Z0) 

CALL ZZLINE( ZO) 

CALL ZZLINE(4) 

WRITE(6,200) N7 

200 FORMAT(///10X,22HSTORAGE NEEDED EQUALS 18) 

IF(N7.GT.AMAX) STOP 333 
CALL ZZLINE(3) 

WRITE(6,300) 

300 FORMAT(//10X,8HVARIABLE,5X,4HMEAN,11X,8HVARIANCE ) 

400 FORMAT(10X,I3,5X,G15.4,G15.4) 

DO 500 I=1,MPM1 
CALL ZZLINE(l) 

WRITE(6,400) I,A(N2+I-1) ,A(N3+I-1) 

500 CONTINUE 
RETURN 
END 

SUBROUTINE PDS ( BUFF , MP , KBUFF , A , XI , XI P 1 , MPMl , DELTT , 

1 XM ,X2M , FREK , LP2 , IWX , L02P 1 , IF02P 1 , IDIV , 

2 RA,RX1,RX1P1,FREKX,L0G) 

DIMENSION BUFF(MP , KBUFF) ,X1 (L02P1 ,MPM1 ) ,X1P1 (L02P1 ,MPM1 ) , 

1 XM( MPMl ) , X2M( MPMl ) , FREK( IF02P 1 , MPMl ) , A( L02P 1 ,MPM1 ) , 

2 RXl (LP2, MPMl ) ,RX1P1(LP2, MPMl) ,RA(LP2, MPMl) ,FREKX(1) 
COMPLEX X1,X1P1,A 

LOGICAL LOG 
C 

C INITIALIZE VARIABLES 

C 

DO 9 1=1, MPMl 
XM(I)=0. 

X2M(I)=0. 

9 CONTINUE 
L=LP2-2 
L02=L02P1-1 
N=0 

IERR=UNIT(8) 
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REWIND 8 
15 CONTINUE 
C 

C CALCULATE MEAN AND VARIANCE 

C 

CALL INOUT(BUFF,MP,KBUFF,l,8,IERR) 

IF(IERR.EQ.O) GO TO 10 
K=LENGTH(8)/MP 
DO 11 J=1,MPM1 
DO 11 1=1, K 

XM( J ) =XM( J )+BUFF ( J+1 , 1) 

X2M(J)=X2M(J)+BUFF(J+1 ,I)**2 

11 CONTINUE 
N=N+K 

GO TO 15 
10 CONTINUE 

DO 12 J=1,MPM1 

XM(J)=XM(J)/N 

X2M(J)=X2M(J)/N-XM(J)**2 

12 CONTINUE 
IERR=UNIT(8) 

REWIND 8 

C 

C K IS THE NUMBER OF RECORDS OF LENGTH L/2 ON FILE 8 

C 

K=N/L02 

N=0 

C 

c LOAD XI AND SET INITIAL A=0 

C 

25 CONTINUE 

DO 30 1=1, LP2 
IF(I.GT.L02) GO TO 35 
IF(N.NE.O) GO TO 40 
CALL INOUT(BUFF,MP,KBUFF,l,8,IERR) 

IF(IERR.EQ.O) GO TO 100 

N=LENGTH(8)/MP 

NI=0 

40 CONTINUE 

NI=NI+1 

DO 31 J=1,MPM1 
RX1(I,J)=BUFF(J+1,NI)-XM(J) 

RA(I,J)=0. 

31 CONTINUE 

IF(NI.EQ.N) N=0 
GO TO 30 
35 CONTINUE 

C 

C SET RX1(I)=0 FOR I L/2 

C 

DO 32 J=1,MPM1 
RX1(I,J)=0. 
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RA(I,J)=0. 

32 CONTINUE 

30 CONTINUE 

DO 36 I=1,MPM1 
CALL FFA(X1(1,I),L) 

36 CONTINUE 

C 

C START ALGORITHM 

C 

DO 50 1=1, K 
IF(I.EQ.K) GO TO 55 
C 

C READ XlPl 

C 

DO 51 11=1 ,L 

IF(II.GT.L02) GO TO 52 

IF(N.NE.O) GO TO 53 

CALL IN0UT(BUFF,MP,KBUFF,1,8,IERR) 

IF(IERR.EQ.O) GO TO 100 

N=LENGTH(8)/MP 

NI=0 

53 CONTINUE 
NI=NI+1 

DO 54 J=1,MPM1 

RX1P1(II,J)=BUFF(J+1,NI)-XM(J) 

54 CONTINUE 
IF(NI.EQ.N) N=0 
GO TO 51 

52 CONTINUE 

DO 61 J=1,MPM1 
RX1P1(II,J)=0. 

61 CONTINUE 

51 CONTINUE 

C 

C CALCULATE FFT 

C 

DO 56 J=1,MPM1 
CALL FFA(X1P1(1,J),L) 

56 CONTINUE 
GO TO 57 

C 

C X1P1(K)=0. 

C 

55 CONTINUE 

DO 58 J=1,MPM1 
DO 58 11=1, LP2 
RX1P1(II,J)=0. 

58 CONTINUE 

57 CONTINUE 
C 

C CALCULATE A 

C 
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DO 60 J=1,MPM1 
DO 60 II=1,L02P1 

A(II,J)=A(II,J)+Xl(II,J)*(CONJG(Xl(II,J))+(-l)**(II-l) 

1 *C0NJG(X1P1(II,J))) 

60 CONTINUE 
C 

C X1=X1P1 AND REPEAT MAJOR LOOP 

C 

DO 50 J=1,MPM1 
DO 50 II=1,L02P1 
X1(II,J)=X1P1(II,J) 

50 CONTINUE 

C 

C SHIFT REAL PART OF A INTO SINGLE REAL ARRAY 

C 

DO 82 J=1,MPM1 
DO 82 I=2,L02P1 
RA(I,J)=REAL(A(I,J)) 

82 CONTINUE 

C 

C TAKE INVERSE FFT OF A 

C 

DO 80 J=1,MPM1 

CALL IFTSYM(RA(1,J) ,L,X1) 

80 CONTINUE 

C 

C DIVIDE A BY SIMGA**2 AND K*L/2 

C RA IS NOW THE ESTIMATE OF THE AUTOCOVARIANCE FUNCTION 

C 

XN=K*L02 

DO 90 K=1,MPM1 

DO 90 I=1,L02P1 

RA( I ,K)=RA( I ,K) / (XN*X2M(K) ) 

90 CONTINUE 
C 

C APPLY WINDOW 

C 

IFREK=(IF02P1-1)*2 

CALL WIND0W(RA,LP2,MPM1,IWX,FREK,IF02P1) 

C 

C TAKE FFT TO GET PDS 

C 

DO 91 J=1,MPM1 

CALL FFTSYM(FREK(1,J) ,IFREK,FREKX) 

91 CONTINUE 
C 

C OUTPUT PDS AND SCALE 

C 

IERR=UNIT(8) 

REWIND 8 

DF=1./(DELTT*IFREK) 

XN=DELTT/(8 .*ATAN( 1 . ) ) 
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F=DF/2 . 

J=0 

DO 70 I=1,IF02P1 
J=J+1 

BUFF(1,J)=F 

IF(LOG) BUFF(1,J)=AL0G10(F) 

F=F+DF 

DO 71 K=1,MPM1 
KP1=K+1 

BUFF(KP1 , J)=SQRT( ABS(FREK(I ,K)*XN*X2M(K) ) ) 

C 

C IF IDIV=N DIVIDE F(N,J) INTO F(I,J) 

C FOR 1=1 TO MPMl, I.NE.N 

C 

IF(IDIV.EQ.K.OR.IDIV.EQ.O) GO TO 71 
BUFF(KP1,J)=BUFF(KP1,J)/BUFF(IDIV,J) 

IF(LOG) BUFF(KP1,J)=20.*ALOG10(BUFF(KP1,J)) 

7 1 CONTINUE 

IF(J.NE.KBUFF) GO TO 70 
J=0 

CALL INOUT( BUFF , MP , KBUFF ,0,8, lERR) 

70 CONTINUE 

IF(J.NE. KBUFF. AND. J.NE.O) CALL INOUT ( BUFF, MP,J, 0,8, lERR) 

ENDFILE 8 

RETURN 

100 CONTINUE 
PRINT 101 

101 FORMAT(10X,20HBUFFER ERROR ) 

RETURN 

END 

SUBROUTINE WINDOW( A, LP2 , MPMl , IWX , FREK , IF02P 1 ) 

DIMENSION A(LP2, MPMl ),FREK(IF02P1, MPMl) 

REAL HANN,HAMM 

BART( I)=FL0AT(L02P1-I+1 ) /XL02P1 
HANN(I)=.5*(1.+C0S(3.14159*FL0AT(I-1)/XL02P1)) 
HAMM(I)=.54+.46*C0S(3.14159*FL0AT(I-1)/XL02P1) 
BLACK(I)=.42+.5*C0S(3.14159*FL0AT(I-1)/XL02P1) 

1 +.08*COS(6.28319*FLOAT(I-l)/XLO2Pl) 

C 

C SET FREK=0 

C 

L02Pl=LP2/2 

XL02P1=FL0AT(L02P1) 

DO 11 J=1,MPM1 
DO 11 I=1,IF02P1 
FREK(I,J)=0. 

11 CONTINUE 

DO 12 I=1,L02P1 
AMP=1. 

IF(IWX.EQ.2) AMP=BART(I) 

IF(IWX.EQ.3) AMP=HANN(I) 

IF(IWX.EQ.4) AMP=HAMM(I) 
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IF(IWX.EQ.5) AMP=BLACK(I) 

DO 12 J=1,MPM1 
FREK(I,J)=A(I,J)*AMP 
12 CONTINUE 
RETURN 
END 

SUBROUTINE INOUT ( BUFF , N , M , IN , UN , lERR) 


C 

C 


c 

BUFF 

INPUT/OUTPUT ARRAY 

c 



c 

IN 

=1, READ DATA 

c 


OTHERWISE, WRITE DATA 

c 



c 

UN 

UNIT NUMBER 

c 



c 

lERR 

NEGATIVE IF NO DATA ERROR 

c 


ZERO IF EOF 


C 

C 

DIMENSION BUFF(512) 

INTEGER UN 
K=N*M 

IF(K.GT.512) STOP 222 
IF(IN.EQ.l) GO TO 1 
C WRITE DATA 

BUFFER OUT (UN,0) (BUFF( 1 ) ,BUFF(K) ) 
GO TO 2 

1 CONTINUE 

C READ DATA 

BUFFER IN (UN,0) ( BUFF( 1 ) ,BUFF(K) ) 

2 IERR=UNIT(UN) 

IF(IERR.GT.O) STOP 111 
RETURN 

END 
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ACSL TRANSLATOR INPUT 


PROGRAM VAN DER POL'S EQUATION 

" PROGRAM CONSTANTS" 

CONSTANT LA =1.0 , TSTOP =50.0 

,PI = 3.14159 

" STATE INITIAL CONDITIONS" 

CONSTANT XIC =-1.5340 , XDIC = .76552 

" DEFINE COMMUNICATIONS INTERVAL" 

CINTERVAL CINT = 0.063 

" DEFINE CONSTANTS FOR PDS" 

ARRAY A( 10000) 

INTEGER AMAX,MP,L02, WINDOW, FREQCY,DIV 
LOGICAL LOG,SKIPP,SKIPT 

CONSTANT AMAX = 10000 , MP =3 

,L02 = 512 , WINDOW = 5 

,FREQCY = 1024 , DIV = 0 

,W =0. , PDSX = 0. 

,PDSXD =0. , LOG = .FALSE. 

,SKIPP = .TRUE. , SKIPT = .FALSE. 

INITIAL 

IF( SKIPT ) GO TO TERM 
END$ "OF INITIAL" 

DYNAMIC 

DERIVATIVE 

" VAN DER POL'S EQUATION" 

XD = INTEG(LA*(1.0-X**2)*XD-X,XDIC) 

X = INTEG(XD,XIC) 

END$ "OF derivative SECTION" 

TERMT(T.GT. TSTOP) 

END$ "OF DYNAMIC SECTION" 

TERMINAL 
TERM.. CONTINUE 

IF( .NOT. SKIPP ) ... 

CALL PDSMAIN( A , MP , CINT , L02 , WINDOW , FREQCY , AMAX , DIV , LOG) 
END$ "OF TERMINAL" 

END$ "OF PROGRAM" 
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ACSL RUN-TIME EXECUTIVE INPUT 


PREPAR T,X,XD 
SET CALPLT=.T. 

START 

PLOT "XHI"=TSTOP , X , "LO"=-5 . 0 , "Hl"=5 . 0 
PLOT "XAXIS”=X,"XLO"=-5.0,"XHI"=5.0, ... 

XD , ”HI''=5 . 0 , "LO"=-5 . 0 
SET SKIPT=.T.,SKIPP=.F. 

START 

PREPAR "CLEAR" ,W,PDSX,PDSXD 
PLOT "XAXIS"=W,"XL0"=0.0,"XHI"=0.72, ... 

PDSX , "HI"=3 .0 , "L0"=0 . 


END 
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